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SUMMARY 

The status of an investigation of four numerical techniques for the time -dependent 
compressible Navier-Stokes equations is presented. Results for free shear layer calcu- 
lations in the Reynolds number range from 10^ to 8.1 x 10^ indicate that a sequential 
alterating -direction implicit (ADI) finite -difference procedure requires longer computing 
times to reach steady state than a low -storage hopscotch finite -difference procedure. A 
finite -element method with cubic approximating functions was found to require excessive 
computer storage and computation times. A fourth method, an alternating -direction cubic 
spline technique which is still being tested, is also described. 

INTRODUCTION 

The quasi -parallel assumption successfully used in boundary -layer -type calculations 
is not applicable for many free mixing flows. The complete Navier-Stokes equations must 
usually be solved for flows which have no single dominant flow direction. This paper pre- 
sents the current status of a detailed investigation of several numerical procedures for 
obtaining steady-state solutions for two-dimensional, high Reynolds number, compressible 
free shear flows using the time -asymptotic approach. In particular, the research has 
been directed toward the solution of mixed subsonic -supersonic flow problems. 

Most published numerical solutions of the compressible viscous time -dependent 
Navier-Stokes equations have been for flows with Reynolds numbers much less than 10^. 
Peyret and Viviand (ref. 1) have summarized these solutions through mid-1973. Taylor 
(ref. 2) also analyzed the literature at that time. Most methods up to the time of these 
surveys used explicit difference schemes. Later, Briley and McDonald (ref. 3) and Baum 
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and Ndefo (ref. 4) published alternating-direction implicit (ADI) calculations. The high 
Reynolds number solutions (Victoria and Steiger (ref. 5), Carter (ref. 6), and MacCormack 
(ref. 7)) were all computed with explicit difference schemes. Recently, additional high 
Reynolds number solutions have appeared, e.g., Holst and Tannehill (ref. 8) and Baldwin 
and MacCormack (ref. 9). These solutions were also computed with explicit methods. 

Although conceptually simpler and more easily coded than implicit methods, explicit 
methods are restricted to small time steps relative to the spatial grid size for numerical 
stability. Consequently, such methods require long computatidri 'times to reach a steady r 
state, especially for flows in which a fine mesh has been used such as in regions of high 
shear. For example, the calculation of a shear layer impinging on a blunt body for a > 
Reynolds number of 10^ by Holst, Tannehill, and Rakich (ref. 10) using the MacCormack 
method requires up to 80 min on a CDC 7600 computer. ' 

The methods under investigation are the following: (1) hopscotch (explicit) finite 
difference, (2) alternating -direction implicit (ADI) finite difference, (3) finite element, 
and (4) implicit cubic spline integration. In addition, some calculations have been made 
with the Du Fort-Frankel procedure. The goal of this study is trie development of an effi- 
cient numerical tool to be used in testing fully two-dimensional turbulence models for a 
wide range of free shear flow applications such as interference heating (shock/shear layer 
impingement), separated flows, jet exhaust noise reduction, combustor design, and tangen- 
tial slot injection. This paper summarizes results of calculations for sample mixing 
problems with Reynolds numbers ranging from 10° to 8.1 x 10 . The procedures are 
compared with respect to their accuracy, computer storage requirements, ease of imple- 
mentation, and total time to steady state for computation of sample problems. 


SYMBOLS 

c speed of sound 

Dj diameter of jet * ■ ■ 

f ,g general functions 

H enthalpy 

L differential operator 

M Mach number 

M^ second derivative of S(x^) 
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p 


density 



respectively 
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B -spline function in equation (14) 


Subscripts: 


i,l 


index denoting grid point spatial location 


nodal, index in finite -element mesh 


stagnation condition 


t,x,y 


derivative with respect to time, x-direction, and y-direction 


Superscripts: 


n,* 


index denoting time level 


A bar over a symbol denotes a dimensional quantity. An arrow over a symbol 
denotes a vector quantity. 

PROBLEM DEFINITION 


To provide a basis for comparison of the numerical procedures, a set of standard 
test problems was selected. 


Sample Problems 

Figure 1(a) shows the mixing problem (case 1) originally chosen for use as the 
standard sample problem. This flow is the mixing of a two-dimensional laminar super- 
sonic (Mach 3) jet and a laminar subsonic flow normal to the jet axis. The peripheral 
velocity Vp is higher in magnitude than the normal velocity component arising from 
natural entrainment of the resultant free shear layer for the same jet issuing into quies- 
cent surroundings. As shown in figure 1(a), the solution domain does not extend to infin- 
ity in either the positive streamwise or normal directions. The peripheral flow is applied 
one or more jet diameters above the corner of the wall, and the calculation is truncated 
one jet diameter downstream from the jet exit plane. This problem thus embodies some 
complicating factors which are often unavoidable in computations of flow fields for real 
vehicles, e.g., a sharp corner and the artificial downstream boundary with a significant 
portion of subsonic outflow. Since the individual effect's of these factors are difficult to 
isolate in the computation of such a flow field, calculations were also made for the related 
problem, mixing of two parallel streams, shown in figure 1(b). The computational region 
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begins downstream of the base of the infinitely thin splitter plates. Such a calculation 
obviously does not require the full Navier -Stokes equations, since solutions can be 
obtained with the usual quasi -parallel approach (boundary -layer equations with free shear 
flow boundary conditions). However, since steady -state calculations could be made with 
another method, solutions were obtained for comparison with computed Navier -Stokes 
results. The mixing of a subsonic stream and a supersonic stream (case 2) shown in 
figure 1(b) was chosen for the study of subsonic' boundary conditions. Calculations were 
also made for the mixing of two supersonic (Mach 3 and Mach 1.68) streams (case 3), a 
flow free from subsonic boundary problems. 

Governing Equations 

The governing equations can be written in nonconservative forms as follows: 
Continuity 


p t + pv y +. pu x + vp y + up x = 0 


( 1 ) 


x -momentum 


pu t ♦ pvu y ♦ pu % = -Px+55^ K)„ - 3 + + v *)] y (2) 

y -momentum 


pv t + pw y + PUV X = -p y + 3j^(l*Vy) y 

These equations are nondimensionalized with 
flow conditions, i.e., 
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The pressure was evaluated by means of the perfect gas equation of state 

■ • ' -7 

P = pRT (4) 

where R = Air was the test gas. Only laminar (molecular) viscous effects were 

considered, the Sutherland law being used to express the viscosity as a function of temper- 
ature 

Id = T 

•t : ■ 

To simplify the system of governing equations and' to reduce required machine storage, 
a constant total temperature of 530° R (294 K) was assumed. Calculations for a Mach 3 
jet into still air with the quasi -parallel code of Oh (ref. 11), which included the energy 
equation, showed that the total enthalpy varied less than- 5 percent throughout the mixing-/ 
region from the constant value assumed in other /calculations. This small variation had* 
a negligible effect on the other flow parameters. As a result of this assumption, the tem- 
perature could be evaluated by the algebraic relationship K 

/ ‘ * ’ .V 1 - ■ • •• • * • 1 . ’ Ui'Ur. 

T = 1 - u 2 - v 2 . . • / (6) 

which eliminated the need for solving the complete energy equation. Constant static . - v 
pressure was assumed in all calculations to generate initial values of density using 
equations (4) and (6) along with the given initial velocities. The linearized version of < 
equations (1) to (6) with the viscous terms neglected has been shown by Gottlieb and 
Gustafsson (ref. 12) to be well -posed for the initial value problem. 

DESCRIPTION OF NUMERICAL PROCEDURES 


3/2 


1 + 


198 . 6/7 


T + 


198. 6/ T £ 


(5) 


Hopscotch 

The hopscotch method is a two-step explicit procedure which was shown to be uncon - 
ditionally stable for the diffusion equation by Gourlay (ref. 13). It was used by Scala 
and Gordon (ref. 14) for compressible viscous calculations of low Reynolds number flow' 

• •>/*,.■'* r* - r i > :*r . ,rff 

around a circular cylinder, and it has been applied to hyperbolic systems with shocks by 
Gourlay and Morris (ref. 15). 

. . Figure 2 shows the pattern of the two sweeps., Consider, for example, the equation 


u t ^xx + u yy 


(7) 
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n+1 


With forward time and centered space differencing, u. j is computed at each time step 
•at the nodes for which i + j + n is even (marked with circles in fig. 2) during the first 
sweep through the mesh with the equation 


= «?. + At 


1,3 


1,3 


n 0 n n 

Vi.i - 2 u i.i tu i-i,i 


(Ax)' 


+ At 


u? 1 . * - 2u? + ii? • i 

i,3+l i,] i,]-l 


(Ay y 


( 8 ) 


This sweep is fully explicit. In the second sweep at this time step, 

n+i' 

"i,3 i-l,j 


: n+l • n ' . . 
Uj j = u- . + At 
. «I: 1,3 1,3. i 


Viz -*&+■** 


(Ax)' 


+. At 


u n+ i, -2u n+1 + u n+ \ 

1,3+1 1,3 1,3-1 


(Ay) 


2 : 


(9) 


at the nodes (marked with squares in fig. 2) for . which i + j+'n is odd. .This sweep is 
implicit in the sense that the values in the computation of the spatial derivatives are at >. 
the new time level n + 1; However, this implicitness does not require the reduction of 
a matrix, since these values were computed during the.first sweep. Differencing which 
does not fit into this pattern, such as a five -point difference for u^ using values of 
U&. and u^j, requires special consideration. The conventional nine -point differ- 
ence analog for cross -derivative terms must receive treatment which usually requires 
the reduction of a matrix. The computational efficiency of the hopscotch procedure is (y 
thus reduced. For the full Navier-Stokes equations, hopscotch has no cell Reynolds num- 

Ay 

ber limitation, but the maximum time step is limited by the condition, At ^ 


For the present application, a sufficient condition for stability is At = 


u + v + 2c 
Ax 


u + v + 


\j2c 


The hopscotch version derived for the present investigation is a low-storage pro- 
cedure (one array per dependent variable). The equations were linearized by lagging the 
nonlinear coefficients. On the second sweep, values at time n + 1 are used only where 
available. This lagging eliminates the need for matrix reductions and thereby simplifies 
the coding, maintains the low storage, and minimizes CPU time^per nodal point. Gottlieb 
and Gustafsson (ref. 12), considering the convective terms only, have analyzed the stability 
of this version of hopscotch with the lagging of some values and have found its stability to 
be identical to that of the original hopscotch method. The method is different, however, 
when the diffusion terms are included. The stability limit which was derived from the 
advection terms is not changed for the range of Reynolds number considered in the present 
investigation. The lagging of values used to compute the viscous terms introduces slightly 
more second-order dissipation than in the original hopscotch method. The new procedure 
is formally not consistent with the time -dependent problem; however, the extra error term 
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introduced is very small for large Reynolds numbers (such as 10^). In the present appli- 
cation the interest is not in the transient but in the steady -state solution; therefore, this 
error term, which goes to zero at steady state, has no detrimental effect. 


During each sweep the x -momentum, y -momentum, and continuity equations are 
solved sequentially at each nodal point with the boundary values then being updated at the 
end of each time step. To illustrate the present version of hopscotch, the differencing of 
the x -momentum equation is as follows: 

First sweep 



+ 2C 


r M n . + M “ \/u n . -u" 
^1,3 *1-1,3 1,3 1-lJ 
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xt n 
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.i 


+ 2F 


r pP . + m” . tWu*. -u n . 
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,n 


+ EM ?i(S) i * +FfjL ii i(i* 

i,J\ax/y 1,3 1 V a ^i > j-l 


(10) 


where A, B, C, D, E, and F are coefficients arising from' the differencing. 


Second sweep 
.n+1 
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= u i,r At< 


u n .m 

1,3 lax 


n+1 

u 





n+1 


(Equation continued on next page) 


444 



3 N Re p ”j 
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Alternating-Direction Implicit Method 

The alternating-direction implicit (ADI) technique developed by Peaceman and 
Rachford (ref. 16) is a two-step procedure requiring reduction of tridiagonal matrices for 
which an efficient solution algorithm, the Thomas algorithm (ref. 17), exists. The method 
was originally applied to the two-dimensional heat conduction equation in reference 16 and 
later to a system of hyperbolic equations by Gourlay and Mitchell (ref. 18). For both of 
these model problems, it was shown to possess unconditional stability. The method, how- 
ever, has not been extensively applied to the compressible Navier -Stokes equations. In 
1966, Polezhaev (ref. 19) obtained solutions for a natural convection problem. His ADI 
method removed the diffusion time -step limitation; however, he found experimentally that 
the time step was still limited to the usual maximum explicit value. In 1973, Baum and 
Ndefo (ref. 4) published a two-dimensional implicit method based on the Peaceman- 
Rachford procedure. The Baum -Ndefo method iteratively solves nonlinear difference 
equations as a sequence of linear equations using a quasi -linearization technique. In a 
one -dimensional calculation of shock structure, the method was found to be stable for 
Courant numbers as large as 10. However, reference 4 does not consider the full Navier- 
Stokes equations. Later in 1973, Briley and McDonald (ref. 3) presented a method based 
on a fully implicit backward time difference scheme in which nonlinearities at the implicit 
time level are linearized by a Taylor’s series expansion about the known time level. The 
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resulting system of multidimensional coupled linear difference equations is solved with a 
noniterative Douglas -Gunn ADI approach. The method was shown to be stable for very 
large Courant numbers in' calculation of three-dimensional subsonic flow in a straight 
duct with rectangular cross section. For a flow with Mach number of 0.044 and a Reynolds 
number of 60, stable solutions were obtained for Courant numbers up to 1250. For a Mach 
number of 0.5 and a Reynolds number of 600, the time step was gradually increased as the 
solution progressed, resulting in an average Courant number of 73i ".Thus, ; the actual 
Courant number decreases with increasing Reynolds number, perhaps because of diagonal 
dominance problems as discussed in reference 3. The computational effort per time step 
was reported to be twice that of; most explicit; methods.,- , . 

In the ADI procedure used in the present investigation, a sequential solution of the ,! 
difference equations is obtained for , each row during the first one -half time step (horizon- 
tal sweep) and for each column during the second one -half time step (vertical sweep). All 
spatial derivatives were approximated by centered finite differences; time derivatives, by 
backward differences. The nonlinear coefficients in the convective terms were lagged one- 
half time step. In addition, the pressure terms and cross -derivative terms were treated 
explicitly in each sweep. The temperature and viscosity were updated, for the entire field 
after each sweep. The order of solution for each row and column is (l) x-momentum equa- 
tion, (2) y-momentum equation, and (3) continuity equation. The solution is then marched 
to steady state without iteration. This ADI formulation requires two storage arrays for 
u, v, p, and p and one for T. 

To illustrate the ADI method, the finite -difference form of the x-momentum equation 
is shown. For the horizontal sweep, from time level n to an intermediate time denoted 
by *, 




, 4 

"fi+ 1,3 + M u) 

l,i 

V 2 A x J 

3N„ Ax 
Re 

\ 2 / 

\ Ax / 


(Equation continued on next page) 
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Finite -Element Method 

The finite -element method has been used extensively for the numerical solution of 
structural mechanics problems for a number of years; however, the procedure has only 
recently been applied to fluid mechanics problems. (See pp. 240-257 of ref. 20.) Using a 
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stream -function — vorticity approach with linear elements, Baker (ref. 21) has developed 
an algorithm for steady viscous compressible flows. Solutions have been obtained for 
Reynolds numbers up to 7750. The present method appears to be the first finite -element 
procedure for the compressible Navier-Stokes 'equations in primitive variables. 


The solution algorithm uses a Galerkin method with leapfrog time integration. The 


solution domain -is discretized with triangular elements with bicubic trial functions. Fig- 
ure 3 shows a finite -element mesh for case 1. As an example of a trial function, the den- 


sity is of the form 
:r '' P(x,y,t) 

. •» i - • 

* i * * , • , * 


10 

‘ ^ Pj(t) '0j(x,y) 

* J=1 ,5 ' *■ •" 




where . J. is the nodal index. and the. ,<pj _ are the so-called B-spline basis functions which 
are piecewise cubic over the problem domain. 


The unknown parameters in each trial function are the flow variable function values 
(p, u, and v) and their first partial derivatives (p x , p y , uy, . v x , and v y J at the 

triangle vertices and the function values alone at the triangle centroid. In the Galerkin 
approach, the weighted residuals formed by using the weights (pj are set equal to zero. 
This yields a set of algebraic equations for the nodal values. Thus, if the governing equa- 
tions are of the form 


L(w) = 0 '/ 

n • ■ * « . 

where w is a general function, the Galerkin approach yields a set of equations 

Jj’ V"> * 0 <15 > 

Solution 

domain 

The time discretization scheme is similar to the Crank -Nicolson Galerkin method 
described by Douglas and Dupont (ref. 22) . Centered time differences over two time 
steps, n - 1 to n + 1, and the averaging of space derivatives over times ' n - 1 and 
n + 1 yield second-order time truncation error. This spatial averaging also elimin- 
ates nonlinearities in the resulting implicit system of difference equations. The system 
of determining equations has the following form: 


Continuity 


n ^ n+1 = FD n,n_1 


D p 


(16) 
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Momentum equations 



(17) 


■ where ^p"' is the unknown density vector/ '{? '• and ! v ' are vectors of unknowns from the 
x- and y -momentum equations, ZZ, ZR, RZ, RR, and D are nonsymmetric matrices 
with varying bandwidth, and FD, FU, and FV are vectors of known quantity. 


The matrices in these equations are assembled at, each time step. The continuity 
equation is solved with a standard triangular decomposition method which takes advantage 
of matrix sparseness. The momentum equations are solved with a unique block iterative 


LU Solver developed during the present investigation. 


. , t • S . , I , ■ 

Despite large matrices and the 


accompanying problem of efficient data management; the use of cubic elements yields 
fourth-order spatial discretization error. Cubic elements also, allow exact incorporation 
of first -derivative boundary conditions, unlike finite -difference methods which require a 
discretization. In addition, the triangular mesh allows the method to be easily adapted to 
nonrectangular solution domains. . 


Cubic Spline Integration Method 

The potential of a cubic spline collocation procedure for the numerical solution of 
partial differential equations has been demonstrated by Rubin and Graves (ref. 23) for 
several model problems. This use of a cubic spline approximation for the evaluation of 
spatial gradients provides a highly efficient and accurate procedure for computation with 
a nonuniform mesh (which is necessary for high Reynolds number calculations in the phys- 
ical plane) and/or curvilinear boundaries. The basic spline approximation leads to a 
second-order accurate expression for second derivatives, e.g., the diffusion terms in the 
momentum equations, for both a uniform mesh and an arbitrarily nonuniform mesh. First 
derivatives, i.e., the convective terms, are third-order accurate with a nonuniform mesh 
and fourth-order accurate with a uniform mesh. With a three -point finite -difference 
approximation, the order of the truncation error is significantly decreased with even a 
moderate variation in the mesh spacing (ref. 24). Thus, the spline procedure is more 
accurate than the usual finite -difference procedures for nonuniform grids. The spline 
method also allows accurate interpolation if grid realinement becomes necessary. 

In addition, first- and second-derivative boundary conditions can be applied more 
accurately and more easily than with conventional finite -difference methods, since dis- 
cretization is unnecessary. Unlike the finite -element or other Galerkin procedures, the 
evaluation of quadratures, which are generally not tridiagonal, is unnecessary. 
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In reference 23, Rubin and Graves present a detailed discussion of the general spline 
formulation and methodology for solving second-order quasi-linear partial differential 
equations. Therefore, only a brief description of the general cubic spline procedure is 
presented in this paper. 

A cubic spline S(x) is a continuous function which has continuous first and second 
derivatives on aminterval ; a <-x < b ‘(a : and’ b ? are two arbitrary points) and corresponds 
to a cubic polynomial in each subinterval Xj_j 5 x = Xj'. ! The mesh 'spacing hv •is defined* 
by \ = x i - x i _ 1 . 

; • The following tridiagonal formulas are obtained by enforcing the continuity require- 
ments at the collocation points xr. 


h, h { + ^ . h, i Ui.i - Uj Ui - Ui , 

8 “»-» * 3 M i + 6 "141 ’ X.1 ‘ • . h, . , . 

(18) 

‘ ♦ *£ + 1 V + „ * + fckli 

h i Vl/ Vi • h^ h 2 

' > . . ' i* 

(19) 

whereat x = Xj, S^x^ = u^ = m A , and The following useful relation- 

ships also exist between the first and second derivatives: 

m i+l - ra i = ^T i ( M i + M i+l) 

% 

(20) 

hi hi U: - uj i 

■ “t = 3 M i + 6 “l-l + h. 

(21) 

m — ^i+1 i\/f ^i+1 A/r i u i+l u i 

m i • 3 M i ‘ 6 M i+1 + h i+1 ' 

(22) 


For a governing partial differential equation of the form 

u t = f ( u > u x> u xx) ( 23 ) 

the approximate solution is found by considering the solution of 



where the time derivative is discretized in the usual finite -difference manner, i.e., ; 

(“t), - ^-sr- 2 - < 25 > 


As : an-exarnple , , consider- an implicit solution of , the linearized Burgers’, equation, which • 
has the general form of the momentum: equations •• 


r * ’ J 


u t+ ,uu x = 


UU 


XX 


I .! ' 


where 



(26) 




u = u(x,t) V = I^(x,t) , ... ' ‘ ! : : 

The approximation for this equation becomes 

. * , • • C . , 

U™ 1 = u!" - 4t(u" 1 m| ,+l ) + 4t(^ +1 Mf +1 ) (27) 

With the spline relations (18) and (19), a system of 3N equations is generated for 
3(N + 2) unknowns. This system can be written as 




( 28 ) 


where V i = Ju^m^My and A, B, C, and D are 3x3 coefficient matrices. Initial 
conditions are prescribed so that u(x,0) = g(x). Equations (20) to (22) can be used, if 
necessary, to relate information at the boundaries and provide a closed system which can 
then be solved by the standard tridiagonal algorithm. 

An alternate procedure can be derived by substituting u- and m. as functions 

has the form 

(i = 1 , . . . , N) (29) 


of Mj. The resulting tridiagonal system for 


a i M ui +b i M r +1 +c i M w= d i 


This procedure is being used in the present application for the two momentum equations. 
If the partial differential equation to be solved has no second -derivative terms (e.g., the 
continuity equation) , a tridiagonal system of equations in terms of m^ can also be found. 
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For the two-dimensional Navier -Stokes equations, a spline ADI procedure has been 
used. This two-step method applies the spline procedure to each of the one -half step ADI 
equations. The cross derivatives are found by using equation (19) with the cross deriva- 
tives being m i and the appropriate first derivatives replacing u^. The three governing 
equations are solved sequentially at each row and column during the horizontal and verti- 
cal sweeps, respectively. 

; ‘ The boundary conditions for u^, u y y,' v^, v y y, p x , and py are found by eval- 

uating the appropriate governing equation at the boundary with the time derivative set equal 
to zero, to give in effect steady-state boundary conditions. The initial values of the second 
derivatives of u and v and first derivatives of p are obtained by fitting cubic splines 
to the given initial function values. 

COMPUTATIONAL RESULTS 

Subsonic Boundary Conditions 

One of the major difficulties associated with the computation of case 1 is proper 
specification of boundary conditions for the region of subsonic flow, i.e., for the subsonic 
portion of the inflow jet profile, the peripheral inflow, and the subsonic portion of the down- 
stream boundary. This boundary -condition problem was therefore studied for case 2 with 
N Re = 8.1 x 10 4 using hopscotch as well as a second-order Du Fort-Frankel procedure 
described by Gottlieb and Gustafsson in reference 25. 

The mathematical analysis of boundary -condition specification by Gottlieb and 
Gustafsson (ref. 12) formed the basis for this study. At the left subsonic inflow boundary 
(see fig. 4) , the analysis indicated that two of the three dependent variables (u, v, and p) 
must be specified. Since v was itself a characteristic variable in the x-direction, it had 
to be one of the two specified functions; u was the logical choice for the second. The 
density boundary condition was chosen to be p x = 0. No difficulties were encountered in 
any of the calculations with this set of inflow boundary conditions. At the upper inflow 
boundary, u is a characteristic variable in the y -direction; therefore, again u and v * ** • 
were specified and p y = 0 was selected as the third boundary condition. This combina- 
tion created no numerical difficulties in any calculations. However; the combination oft^*' r 
p and u specified with v y = 0 usually led to erroneous values for v, especially in the 
region near the upper boundary where positive values of v, indicating outflow, occurred. 

At the subsonic outflow boundary the one -dimensional analysis indicated that one 
function value, either p or u, must be specified. Of course such a boundary condition 
is not convenient for most applications since downstream function values are generally 
not known a priori. Figure 5 shows the results of calculations using hopscotch with three 
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different subsonic downstream boundary conditions. The boundary conditions used in the 
subsonic region are indicated in figure 4. 

The initial flow field was obtained by setting the downstream (outflow) boundary 
values for u and v equal to one -half of their steady -state values (obtained from the 
parabolic code described in ref. 11) and linearly interpolating to obtain values at interior 
nodes. Computed steady-state profiles of the streamwise and normal velocity components 
u and v, respectively, at the downstream boundary are compared with results obtained 
with the parabolic technique. The streamwise component is accurately predicted for all 
three boundary conditions; however, only the specification of p 'gives a smooth and accu- 
rate v,. profile. Specifying u produces large oscillations in the v ; profile in the vis- 
cous region. These oscillations may be critical in turbulent flows when the turbulence 
model is locally a function of 9v/ 9y. The least accurate results are obtained for linear 
extrapolation of all three function values. 

The results of calculations with the Du Fort-Frankel procedure (see ref. 12) were 
identical to the hopscotch results with the exception that the linear extrapolation had to be 
altered to obtain converged solutions. Extrapolation of values at time level n + 1 to 
obtain boundary values does not work for any degree of extrapolation (linear, quadratic, 
etc.). Linear extrapolation of the form 

f? +1 i = 2ff li-fj 1 ' 1 _ 2i (30) 

^ax” max a >j max 


for p, u, and v, where i max is the outflow boundary, gave results which converged to 
the correct steady state. The boundary condition 


jn+1 _ fn 

^nax’^ ^ax" 1 ^ 


(31) 


which has been shown to be stable for scalar hyperbolic equations for the pure leapfrog 
scheme in reference 26, also gave good results. Using both values at n + 1 in equa- 
tion (31) results in an unstable condition. ; 

Parallel Mixing Calculations 

Calculations of case 3 were made with the hopscotch and ADI methods for Np e =10^ 
and 5 x 10^. The supersonic inflow, supersonic outflow, and upper inflow boundary condi- 
tions shown in figure 4 were used. The computed steady state u and v velocity pro- 
files and pressure profiles at x = 0.15 for N Re = 10^ are compared in figure 6 with the 
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corresponding calculation using a parabolic method. The ADI and hopscotch results are 
virtually identical with both procedures accurately predicting u and v. For this grid ' 
spacing (Ax = Ay = 0.025), the pressure shows high-frequency oscillations in the viscous : 
mixing region, although the maximum oscillation is only about 2 percent of the correct 
value, ., For the initial flow field the y = 0 inflow profiles were also specified at all other 
Xrstations downstream. The time. step in the hopscotch calculation was 0.9 of the maxi- 

... , . , > . i i-ir ’ t-i \ i /V. , iS • > - : • . <> - ; Ujrt, f ' 

mum allowed for stability, i.e., the Courant number was 0.9, Identical solutions 

" " ‘ ‘ " ’ • * ‘ 1 '- I 1 * • . ■ {«»..'■: I ' - < . . ! . i ^O. . . , . , :i ' i' ;■ i > h fi'tc.j I'.r ' 

were also obtained with ADI for Nq q = 6 but with no decrease in the total number of step's 

to steady state. The solution diverged for N^, Q =12. 

iff’ » . , , • , , \ 

‘ At both Nq 0 = 6 and 12, the large time step. was used for the entire solution. The 
probable cause of divergence at N Co = 12 is roundoff error from the tridiagonal matrix 


inversion occurring; when, the coefficient matrix did not possess, diagonal dominance , a 
sufficient. but not, necessary condition for convergence of the matrix reduction. In this 
instance the continuity equation was not diagonally dominant for any row in the horizontal' 

' *' - - .1 * ■■ >■ ; ; ~ > i • i t r ' j 

sweep , and in addition, the two momentum equations lacked diagonal dominance for many 

' -• • •’ • . • • (•.. i ! i ; c. •• , •. 

rows. , , , . . , . ...... ... , ■ ’ ' 

. - • , j j . * •• . . *• •. .. 

Figure T. coiitains the steady-state results at ;x = 0.15 for Nj^ e = 5.0 x 10^. .Fig- 
ure 7(a) shows that again the hopscotch and ADI results are virtually, identical for .u . and ; 
v. With the same grid spacing as in the.previous calculation, u is accurately predicted,! 
whereas the < v profile exhibits an oscillation in. the viscous region. Halving the grid so, 
that ‘ Ax = Ay =• 0.0125 eliminated this oscillation; The pressure profile shown in fig- ... 
ure 7(b) has very small oscillations for both grids with the hopscotch pressure varying less 
from the constant value given by the parabolic code than the ADI pressure. For = 0.8, 
with either grid, the ADI method required approximately 5 times as many steps for conver- 
gence to steady state as hopscotch. The solution was considered to be converged when 

f n+l jn-* : 

•=••• : - i l 3 - At ' M = °- Q1 g y ••• >'•'* /. \ ‘ , (32) 

for f = p, u, and v at every, point in the field and g = p, u, and u, respectively! 
Converged ADI solutions were obtained for N^ 0 = 6, but the solution again diverged for 
N^ q =12. Hopscotch solutions were not attempted with Courant numbers significantly 
greater than one. 


•> • Figure 8 shows results of hopscotch calculations for case 3 with N^ e = 8.1 x 10*. 
Profiles of v and p are shown at the downstream boundary x = 0.45. The boundary 
conditions shown in figure 4 were used with function values obtained from the parabolic 
procedure providing the necessary specifications of v at the upper boundary and p at 
the downstream boundary. The initial field was obtained by using steady -state values of 
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all functions (obtained from parabolic code) at the downstream boundary and calculating 
the interior values by linear interpolation. This gives an initial flow field which is a good 
approximation of the steady-state field. For a grid spacing of Ax = Ay = 0.025, the hop- 
scotch solution converged in 1637 time steps. The u velocity component was accurately 
predicted, and figure 8 (a = 0 curve) shows that v also agreed with the parabolic results; 
however, small oscillations occurred in the pressure. It was found that the' pressure could 

■Zi( ill 1 j : . {}}«’•! t ? •> *•>«' < ' ■ if (• t :■ ■ 1 • , i. . . - ,, 

be smoothed by explicitly adding artificial diffusion to the' continuity equation, that is; '• £ ’ J :' i 


iP t ;+,(pu) i: .+ (py)y = a(Ax)fp xx + Q!(Ay) 2 p. 


yy 


i -• . 


. . .* - * 


( 33 ) 




Victoria and Widhopf (ref 1 ’27) also found ’it’ necessary to add artificial diffusion to the'* £ ' • 
continuity equation. For this case, a = 0.1' smoothed the pressure without altering u ' 
or v; moreover, the solution converged in 1150 time steps. (The coefficient of p^ 'is ' 
then approximately 5 to i0 times as large as the average coefficients of the viscous terms 




: • i . 


in the two momentum equations.) Although the solution converged in 1837 time steps for 
a = 1.0, both p and ; v show oscillations. ^ For ; a = 0, a significantly different initial 
flow field was generated by halving u and : v' j . at the downstream edge and then linearly' - 
interpolating for interior values. ! Steady state was reached in 3050 steps with u, v, 
and p found to be identical to the previous results. For this gridj no converged ADI - - - 
results were obtained with or without artificial diffusion. Oscillations in v in the mixing 
region grew with time, and the solution diverged. 


The Du Fort-Frankel procedure with the downstream boundary condition given by 
equation (31) for u and v was compared with Hopscotch which used linear extrapolation 
when the initial flow field was obtained by setting the outflow values of u and v to one- 
half of their steady -state values and linearly interpolating for interior values. Both solu- 
tions converged to steady state at approximately the same nondimensional time, although 
slightly different time steps were used in each method. The maximum difference in the 
u profiles was approximately 3 percent and occurred in the viscous region. This slight 

difference is attributed to the difference in explicit artificial diffusion in the two methods. 

v*’J r " ■ ’ • -- 

Second-order diffusion with a = 0.1 was used in hopscotch, and fourth-order diffusion 
was added to the Du Fort-Frankel procedure. 


Results for Case 1 

Computations were attempted for the original test problem, case 1, with a Reynolds 
number of 8.1 x 10^ with the two finite -difference methods and with the finite -element 
method. A nonuniform grid was used for the finite -difference methods. From the sharp- 
corner, Ay was increased by a factor of 1.05 for each successive spacing in the positive 
and negative y-directions from the smallest value, Ay = 0.005. Thus, grid points were 
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concentrated in the viscous mixing region. The grid also increased in the x-direction by 
a factor of 1.05 (ax^/Ax^ = 1.05) for each' .successive spabing. 

Computations with hopscotch were made with the boundary conditions of figure 4 
along with a no-slip wall condition and constant density assumption near the wall. The 
upper boundary .was approximately six jet diameters above the center line. Figure 9 
shows steady -state hopscotch results for an interior x-location using linear extrapolation 
for density in the subsonic portion of the outflow boundary (extrapolation used since p 
is unknown). The magnitude of the peripheral velocity at the upper boundary was 0.07, 
which is several times greater than the natural entrainment for the same jet issuing into 


still, air.. Since the.upper bqundary effectively models a porous wall with nonuniform mass 
injestion.into the boundary layer, it was necessary to modify the upper boundary conditions 

i * - *»* f -l'/?.*. ’» ‘.Tiw . ,~.l i i ’ y i *.*./;• i*. ^ ' • *“ 

by. setting \i v = 0. The results shown in figure. 9 were obtained with a ='10 after 3 hr 
of computing time on a CDC 6600 computer. As expected from the parallel mixing results, 
the u component is smooth and appears to be qualitatively correct. (There are no known 
experimental data for such a flow with which to compare the computed r e suits. ) : The v 
profile shows the oscillation characteristic 'of using linear extrapolation for the subsonic 
outflow density. For some engineering applications, however, these results may be suffic- 


ient. The local increase in the pressure profile indicates the presence of a weak shock. . 

Fully converged ADI results were not obtained for case 1. With linear density 
extrapolation and the same nonuniform grid- the solutions appeared to be nearly converged 
after approximately 3 hr’ of CPU time on the CDC 6600, but the computations were not 
continued further since the results did not appear- to be better than hopscotch. As with 
hopscotch, the u profiles were smooth and apparently qualitatively correct, although v 
again exhibited spatial oscillations. 

At the present time, converged results have not been obtained with the finite -element 
method for this problem. The major difficulties appear to be the lack of sufficient spatial 
resolution in the viscous region and incorporation of the second -derivative downstream 
continuation boundary conditions. The 103 -triangle mesh currently in use (shown in fig. 3) 


requires excessive machine storage and prohibits a significant increase in resolution. 


The cubic spline algorithm has been coded, but presently no steady -state results 
have been obtained. 


CODE COMPARISON 

For case 1 the hopscotch code requires machine storage of approximately 50000g 
for 3045 node points, whereas the ADI method requires approximately 150000g for the 
same grid. The 103-triangle finite -element mesh which has 301 nodes requires 330000g. 
The cubic spline algorithm will presumably require fewer grid points for accuracy com- 
parable to the finite -difference results. 
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For each At step the CPU time on the CDC 6600 for hopscotch is 1.08 x 10'^ 
sec/node and for the ADI method, 3.74 x 10^ sec/node. The finite -element code requires 
1.93 x 10“1 sec/node to assemble the matrices and even with the fast block iterative 
solver, requires 4.25 x 10"^ sec/node for equation solution at each time step. No data 
are available for the spline code. (The Du Fort-Frankel code requires 1.28 x 10'^ 
sec/no^e for each time step.) / “ ■ " '' i ' i 


CONCLUDING REMARKS 

' ’ * ‘ • ' • f - - • • • i ’■ ’ * 


A study of mixed supersonic -subsonic free shear fiows has shown that correct cal- 
culation of the normal velocity component required specification of the density in' the sub- 
sonic portion of the outflow boundary! The streamwise velocity was, however, 'correctly 
computed even when linear density extrapolation was used for this boundary - ' ' ' * - ■ ' 


For high Reynolds number flows (flows with small viscous terms in the. momentum 
equations), it was necessary to add artificial diffusion in the continuity equation to elimin- 
•. ate oscillations in the static pressure. The .addition of . too much artificial viscosity had 
adverse effects on both the pressure and the. normal component of velocity. 


For the problems considered, the maximum allowable time step for the sequential 
alternating-direction implicit (ADI) procedure was less than 10 times the maximum 
explicit time step. This increase in time step, however, did not significantly improve 
the convergence rate. The hopscotch procedure, with a time step no greater than the 
maximum explicit time step, still converged faster than any of the ADI solutions. A fully 
coupled ADI procedure may allow larger time steps; however, the effect of large steps on 
convergence rate must be investigated. 

The finite -element method with cubic elements appears to have excessive storage 
requirements and computing times. Therefore, as currently formulated, it does not appear 
to be a competitive procedure for high Reynolds number calculations in aerospace vehicle 
analysis. 
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(a) Mixing of laminar supersonic jet with imposed peripheral flow 
normal to jet center line axis. 



(b) Mixing of two parallel flows . 


Figure 1.- Standard sample problems. 


Case 2 


Case 3 
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At time level n 
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Figure 2.- Hopscotch grid. 

Inflow BC 

Wall. 
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Inflow 

BC 


Symmetry BC 

Figure 3.- Solution domain showing finite -element mesh for case 1 with 103 triangles 

with boundary conditions (BC) indicated. 
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Subsonic inflow 



Figure 4.- Schematic of computational domain with best boundary conditions for case 2. 


O Parabolic solution 

p specified for subsonic outflow 

Pxx^ = 0 for subsonic outflow 

u specified for subsonic outflow 
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Figure 5.- Boundary -condition study for subsonic boundary. Np e =\8.1 x 10^; x = 0.45. 
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(a) Velocity profiles; 

Figure 7.- Subsonic -supersonic mixing for N Rp = 5 x 10 3 . x = 0.15; N = 0 . 9 . 



P t + (pu) x + (pv) y = (Ax) 2 a Pxx + (Ay) 2 a P yy 



Figure 8.- Effect of artificial diffusion in continuity equation. 
N Re = 8:1 x 10 4 ; x = 0.45. 



